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SECTION  I 
•INTRODUCTION 

The  MASTER  program  file  contains  programs  and  cross-section 
libraries  used  in  solving  radiation  transport  problems.  This 
program  file  includes  as  subsets  the  FASTER,  BETA,  and  TEMPER 
programs  (refs.  1-8)  and  other  analysis  procedures  which  all 
operate  from  a common  data  base.  The  program  is  organized  to 
take  advantage  of  the  similarities  in  different  transport 
problems.  This  organization  includes  input  processors  which 
accept  general  descriptions  of  geometry,  source  distribution, 
detector  distribution,  and  particle  cross  sections. 

The  program  data  base  features  variable  format  data  input 
with  variable  dimension  data  storage.  The  inpuL  for  data  base 
preparation,  and  for  invoking  analysis  procedures,  are  given  in 
the  Users  Guide.  This  guide  is  maintained  as  an  updatable  card 
image  file  as  part  of  the  program  file.  Control  cards  for 
listing  this  guide,  and  for  executing  a problem,  are  svaranarized 
in  the  Appendix. 

Analysis  capabilities  range  from  simple,  one-dimensional 
problems  through  complex  geometry,  nonlinear,  time-dependent 
problems.  Particle  physics  models  include  primary  and  second- 
ary photons,  electrons,  neutrons,  and  heavy  charged  particles. 
Numerical  methods  include  approximate  transport  kernels,  direct 
numerical  integration,  and  Monte  Carlo  quadratures. 

The  general  analysis  procedures  use  Monte  Carlo  methods. 

For  the  more  simple  analyses,  the  Monte  Carlo  method  can  be 
understood  in  the  context  of  particle  simulation.  These 
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general  analysis  procedures  are  implemented  in  the  program 
structure  by  serial  execution  of  three  modules 

a.  Data  base  preparation;  input  and  interpretation  of 
data  describing  the  transport  problem  (subprogram 
DEFINE) 

b.  Particle  tracking;  simulation  of  particle  interactions 
while  accvimulating  contributions  to  detector  response 
using  physics  models  appropriate  to  the  analysis 
(subprogram  SOLVIT) 

c.  Output;  reduction  and  printout  of  sample  statistics 
for  the  accumulated  detector  response  (subprogram 
MSWER) 

Analysis  procedures  are  described  by  starting  with  general 
models  and  then  proceeding  through  specific  models.  The  inter- 
play of  physics  models  is  described  first.  This  description 
is  given  relative  to  an  electron-bremsstrahlung  problem.  Specific 
analysis  methods  are  then  described  and  their  requirements 
relative  to  the  data  base  established.  Finally,  methods  and 
models  implied  by  the  data  base  are  described. 

Recommendations  for  simplifying  program  use,  for  minimizing 
system  impact  (reducing  computer  core  and  time  requirements), 
and  for  improving  physics  models  are  given  in  the  final  section. 
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The  coordinates  of  particles  are  denoted  by  the  phase  space 
vector  P which  has  the  natural  partitions 


P = x,u,E,t 


(1) 


where  x denotes  position,  u denotes  direction,  E denotes  kinetic 
energy,  and  t denotes  time. 

A transport  problem  is  described  by; 

a.  The  fixed  differential  source,  S*^(P) 

b.  The  fixed  differential  response,  D*^(Q) 

c.  The  microscopic  total  and  differential  cross  sections 
for  interaction  and  the  spatial  distribution  of  material 
mixtures 

The  distribution  of  material  mixtures  will  be  assumed  constant 
over  subsets  of  space,  called  regions,  with  discontinuous  changes 
at  region  boundaries,  called  svirfaces. 

Macroscopic  cross  sections,  both  total  and  differential,  are 
formed  by  n\imber  density  weighting  of  the  microscopic  cross 
sections.  This  weighting  procedure  is  implied  in  the  following 
discussions  where  microscopic  cross-section  notation  is  used. 
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The  effects  of  particle  interactions  are  expressed  in  terms 
of  cross  sections  as  the  no-interaction  kernel  (the  differential 
flux  per  imit  source) 


N(P,Q)  = exp(-(?'r)  (5„(Q  - U(P,r)) 


(2) 


where  O’r  denotes  optical  thickness  (mean-free-paths) 


= /0'(u(p, 

Jo 


r'  ) ) dr' 


(3) 


and  the  single-interaction  kernel  (the  differential  source  per 
tinit  flux) 


I(Q,P)  = (Q,p)  (5p(P  - V(Q,p)) 

dp 


(1+) 


In  general,  (5^(»»»)  is  the  Dirac  delta  function  in  all  compo- 
nents of  the  argument  except  t.  Thus,  U(P,r)  denotes  the  phase 
space  coordinates  of  a particle  which  started  at  P and  has 
traveled  the  distance  r,  and  V(Q,p)  denotes  the  phase  space 
coordinates  of  a particle  which  enters  an  interaction  with 
coordinates  Q and  whose  phase  space  vector  component  q has 
been  changed  to  p by  the  interaction. 
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It  is  siifficient  to  represent  the  solution  of  the  transport 
equation  by  a Neumann  series 


= ys^(P)  N(P,Q)  dP 
S^‘^''(P)  = I(Q,P)  dQ 


D°(Q)  dQ 


where  k = 0,  1,  ...  denotes  the  number  of  particle  interactions, 

denotes  the  integral  response  of  the  detector  from  exactly 
k interactions,  and  the  iteration  process  is  terminated  when 
the  response  from  all  future  interactions  is  deemed  negligible. 

The  Monte  Carlo  method  is  used  to  evaluate  the  Neumann 
series  integrals,  i.e. 


s“  = s''(p^)  / p*(P^) 


= s'*  N(P^,Q) 


/ P^(Qn) 


s‘‘'^’(P)  = I(Q^,P) 


The  discrete  phase  space  vectors  and  are  selected  at 
random  from  the  probability  density  f-unctions  P^(P)  and  P^(Q) 
which  are  defined  such  that 


pJ(P)  0 ^ if  S^(P)  =»  0 


pJ^(Q)  =*  0 if  » 0 


and  normalized  such  that 


= dP  = Jp^'(Q)  dQ 


If  the  probability  density  functions  are  defined  as 


pJ^(P)  = S^(P)  / Js^{?)  dP 


pJ(Q)  = 


/ J 25^(Q)  dQ 


then  the  Monte  Carlo  integration  reduces  to  a direct  simulation 
of  particle  transport  with  non  ab  sorp  t ion  and  none  scape  weighting. 
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Random  Sampling 


The  random  selection  of  discrete  phase  space  vectors  P, 

A 

and  implies  random  selection  of  individual  components  of  these 
vectors  from  marginal/conditional  distributions.  Random  sampling 
techniques  used  for  each  phase  space  vector  component  all  utilize 
random  numbers,  denoted  by  r,  uniformly  distributed  on  the  open 
interval  (0,1). 

Solution  of  inverse  functions  is  used  for  sampling  simple 
continuous  distributions 


>''(x)  = r 

-'-o 


p" (x' ) dx 


’ = I dr’  = 

Jo 


X = P (r) 


For  example,  a point  imiformly  distributed  on  the  open  interval 
(a,b)  is  obtained  as 


dx’  = X - a = r 


b - a 


b - a 


X = a + r (b  - a) 


Table  lookup  of  tabulated  cumulative  distribution  functions 
is  used  for  complicated  distributions 


J 

-'-oo 


(x')  dx*  = 


(pI,  - pt  ) 

^i+1  - ^i 


X = x^  + i-H 


1 (r  - pr  ) 


where  x^  , i = 1,  2,  are  points  at  which  the  cumulative 

distribution  is  tabulated  and  the  interval,  i,  containing  the 
point  X is  determined  as 


p'i  = r ^ p*(x' 

—aa 


) dx’ 


Rejection  techniques  are  used  to  sample  some  analytic  functions. 


i*(x)  = 


10 


where 


■/, 


oo 

g(x)  h(x)  dx 

oo 


h(x)  rs  1 


. /:, 

-Loo 


g(x)  dx 


r 

G(x)  = / g(x* ) dx’ 

J-oo 


(23) 


The  distribution  is  sampled  as 


X = G"^  (r)  if  h(x)  a-  r’ 


(214.) 


where  G"'*(r)  is  the  function  inverse  to  G(x)  and  r and  r'  are 
two  different  random  numbers  on  (0,1).  The  procedure  is  repeated 
iintil  an  x satisfying  the  acceptance  condition  is  selected.  The 
normalizing  constant  C is  the  efficiency  of  the  technique,  i.e., 
the  average  number  of  repetitions  to  obtain  an  acceptable  x. 

Random  sampling  of  summations  is  a particular  interpretation 
of  sampling  continuous  distributions,  e.g.,  discrete  interaction 


events  can  be  selected  at  random  using 


KQ^.P)  = ^ 


P)  \ 

Pu(e) 


e Pk^e) 


e,;  e pj(e)  (P5) 


where  equality  is  in  the  expected  (average)  value  sense,  the 
subscript  e denotes  different  events,  and  the  probability 
density  function  has  the  properties 


if 


/ 


dP 


(26) 


The  particular  event  e^^  is  selected  using  a random  number  r on 
(0,1  ) as 


Sample  Problem 


To  illustrate  the  Neumann  series  integration,  a simple 
geometry  is  assumed.  It  will  suffice  to  consider  a semi-infinite 
plate  containing  a single  material# 


The  description  of  the  problem  geometry  requires  the  selection 
of  a reference  coordinate  frame.  For  convenience,  the  material 
slab  is  oriented  perpendicular  to  the  z-axis  with  its  bottom  on 
the  x-y  plane  such  that  its  thickness  is  measured  in  the  +z 
direction.  Letting  h denote  the  slab  thickness,  the  material 
is  bounded  by  two  plane  surfaces 


The  position  of  particles  within  the  material  is  denoted 
generally  by  x and,  more  specific,  by  the  rectangular  coordinates 


Much  more  general  problem  geometries  can  be  described, 
concepts  used  in  more  general  problems  are  discussed  under 
data  base  description. 


A source  of  radiation  is  needed  to  make  the  problem  meaningful. 
A monoenergetic  electron  beam  is  assumed  to  impinge  on  the  bottom 
of  the  slab.  For  simplicity,  the  beam  is  assumed  to  have  no 
breadth  and  to  enter  the  material  slab  at  the  origin.  Thus  the 
coordinates  of  electrons,  as  they  enter  the  slab,  are 


^ = (0,0,0) 

The  initial  direction  of  the  electron  beam  is  given  by  the 
polar  angle  (j)  measured  from  the  z-axis,  and  the  azimuthal 
angle  Q measured  from  the  x-axis.  For  the  described  problem, 
the  initial  azimuth  is  not  critical.  • It  will  come  into  play 
as  the  particle  traverses  the  slab. 

It  is  convenient  to  denote  directions  by  direction  cosines 
relative  to  the  three  coordinate  axes,  i.e.. 


u = (C?  ,/},/)  = (c^,C2,c^)  = CXi  + (29) 

where 


c^  = cos  Q sin  (j) 
C2  = sin  B sin  (j) 
c^  = cos  0 


(30) 


1li 


The  cosine  of  the  polar  angle  will  almost  always  be  used 
rather  than  the  polar  angle. 

In  a more  general  problem,  the  radiation  source  may  have  a 
spatial  extent  in  one  or  more  of  the  three  rectangular  coordi- 
nates, or,  equivalently,  in  one  or  more  of  the  three  spatial 
variables  in  cylindrical  or  spherical  geometry  source  volumes. 

The  radiation  source  may  also  start  in  various  directions 
ranging  from  the  monodirectional  beam  described  here  to 
isotropic  and  even  complicated  distributions  in  the  two  angular 
variables . 

Assuming  the  electron  besim  for  the  sample  problem  is-  mono,-- 

' p 

energetic,  the  initial  particle  energy  is  denoted  by  Eq.  Particle 
energies  are  measured  in  MeV  and  denote  the  kinetic  energy  only. 
This  electron  can  be  visualized  as  a single  particle  or  it  can 
be  said  to  represent  the  entire  beam,  depending  upon  a prefer- 
ence for  relative  or  absolute  calculated  output  units.  The 
m;unber  of  electrons  represented  by  this  one  particle  is  called 
the  weight  and  is  denoted  by  W or  S^. 

For  more  complicated  sources,  a spectrum  of  energies  may  be 
present  in  the  source.  This  spectrum  can  be  described  by  a 
tabulated  distribution  and  would  be  sampled  to  obtain  a single 
mono  energetic  particle  using  a tabulated  distribution  sampling 
technique,  i.e.,  table  look-up. 

Assuming  the  beam  of  particles  enters  the  plate  at  a refer- 
ence time  t = 0,  completes  the  specification  of  the  initial  phase 
space  coordinates. 


Electron  Tracking 
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The  first  step  in  tracking  a particle  is  determining  the 
maximvmi  distance,  and  the  materials  traversed,  for  the  problem 
geometry.  The  simple  slab  geometry  yields  the  general  distance 


s = (h  - z^)  / c^ 


if  c^  =»*  0 


s 


if  c^  0 (31) 


The  next  step  in  particle  tracking  is  determining  whether 
the  particle  has  an  interaction  within  the  defined  geometry 
or  leaks  out  through  one  of  the' bovindaries  of  the  problem  geometry. 

Electrons  have  a large  number  of  interactions  with  the  orbital 
electrons  and  with  the  field  of  the  nucleus  as  they  traverse 
material.  It  is  computationally  inefficient  to  treat  each  of 
these  interactions  explicitly.  Instead,  the  electron  is  allowed 
to  traverse  a small  distance  ignoring  interactions.  At  the  end 
of  the  free  flight,  the  electron  coordinates  are  corrected  for 
the  net  effect  of  the  interactions  ignored  on  the  free  flight. 

This  procedure  is  called  the  condensed  history  method. 

The  first  material  effect  on  electrons  (and  other  charged 
particles)  is  continuous  energy  loss  along  the  free  flight  path. 

The  processes  causing  the  energy  loss  are  inelastic  interactions 
with  the  atomic  electrons  (with  creation  of  secondary  electrons) 
and  radiative  interactions  with  the  field  of  the  nucleus  (with 
creation  of  bremsstrahlung  photons). 
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The  effect  of  these  energy  loss  terms  is  denoted  by  L(E) 
and  for  electrons  and  other  charged  particles  is  called  the 
stopping  power.  Stopping  powers  have  units  of  MeV/cm  and  rep- 
resent the  average  energy  loss  per  unit  distance  traveled  by 
the  particle.  There  is  an  equivalent  quantity  for  neutral  par- 
ticles which  has  the  same  definition  but  is  usually  called  the 
linear  energy  transfer  coefficient  and,  for  photons,  the  energy 
absorption  coefficient. 

In  a general  sense,  L(E)  can  be  calculated  as 


L(E) 


— (E,E>)  (E 
dE' 


E')  dE' 


(32) 


where  d CX/dE'  is  the  probability  (cross  section)  per  unit  path 
length  for  scattering  the  particle  from  energy  E to  E',  per  unit 
scattered  energy,  and  (E  - E' ) is  the  energy  lost  in  the  scat- 
tering process.  Written  in  this  form  there  is  an  implied  svim- 
mation  over  different  reactions  for  a given  element  of  the 
material  and,  for  multiple  element  materials,  an  atom  density 
weighted  sximmation  over  the  elements. 

The  electron  free-flight  path  length  is  determined  using 
the  stopping  power.  The  path  length  is  fixed  by  allowing  the 
electron  to  lose  a fraction,  f,  of  its  current  energy  which 
yields  the  implicit  definition  of  the  free-flight  distance,  d 

f dE» 

d = / (33) 

^ L(E'(s')) 
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After  the  electron  traverses  the  distance  d,  particle  coordi- 
nate updates  are  made 


= E - f E 


^ “ ^-1 


+ d 


where  the  index  k denotes  the  number  of  electron  steps  since 
the  history  was  started. 


From  the  updated  position,  it  can  be  determined  if  the 
particle  escapes  the  material  slab.  If  x^  is  less  than  zero, 

the  electron  has  escaped  from  the  incident  side  and  is  called 
a reflection.  This  cannot  occur  on  the  first  step  of  course. 
If  x^  is  greater  than  the  slab  thickness,  the  particle  has 
escaped  from  the  top  of  the  slab  and  is  said  to  be  transmitted. 

The  time  of  the  particle  is  also  updated  as 


^k  = ^k 


I ds' 

J.  v(E«(s» 


where  v(E* ) is  the  instantaneous  particle  speed  (cm/sec) 


v(E)  = 


E (E  + 2 m 

o ' 


E + m c 
o 


where  c is  the  speed  of  light  (not  to  be  confused  with  any  of 

2 

the  direction  cosines),  m is  the  particle  rest  mass,  and  m c 

o o 

is  the  particle  rest  energy. 


If  the  particle  has  not  escaped  from  the  material  slab, 
the  next  procedure  is  to  account  for  the  interactions  neglected 
on  the  free-flight  path.  These  neglected  interactions  include 
the  net  effect  of  multiple  angular  deflections  and  the  net 
effect  of  energy  losses  less  than  and  greater  than  the  mean  loss. 

The  effect  of  the  angular  deflections  is  measured  relative 
to  the  direction  u,  . The  deflection  is  assumed  equiprobable 

in  azimuth  around  this  fixed  direction  and  varies  with  the  cosine 
of  a polar  angle  measured  relative  to  this  fixed  direction  as 
dictated  by  the  Goudsmit-Saimderson  distribution 

p(/Z)  = ^ A„(d)  p^(/i)  (37) 

n = 0 


An(d)  = 


^ ^ exp  -f  O'n(s')  ds' 

2 L Jo 

cr„(e)  = J [l  - P^,(/Z)] 


(36) 


(39) 


d/i 


where  d Cf/dli  is  the  macroscopic  angular  deflection  crcoS  section 
(per  -unit  deflection  angle  cosine)  and  fJi)  is  the  nth  Legendre 
polynomial. 
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Since  the  angular  deflection  is  sharply  peaked  in  the  forward 
direction,  the  multiple  deflection  distribution  has  a similar, 
but  not  as  severe,  forward  peaked  distribution.  The  cosine  of 
the  deflection  polar  angle,  relative  to  Uj^  , is  obtained  by 

random  sampling  of  a cumulative  distribution  function  tabulated 
on  a fixed  mesh  of  electron  energies  and  deflection  cosines. 

With  the  random  selection  of  an  azimuth  around  , the 

direction  cosines  of  the  electron  in  a coordinate  system  with 
z-axis  along  u^^  are  denoted  by  (cjj , c^,  c!^) . These  direction 

cosines  are  rotated  to  the  reference  coordinate  system  as 


c . = 
J 


3 

z 

i = 1 


°i  ^i,j 


j = 1,  2,  3 (i^O) 


where  the  rotation  matrix  R is  determined  by  the  polar  and 
azimuthal  angles  of  u^  measured  in  the  reference  system 


R = 


cos  cos  ^ 


- sin 


cos  0,  sin 


sin  (7,  cos 


- sin 


sin  0,  sin  cos 


The  same  rotation  procedure  is  used  on  all  particle  reactions, 
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If  bremsstrahlung  photons  are  being  created  in  the  material 
slab,  they  can  be  modeled  as  created  at  random  anywhere  along 
the  free-flight  path  or  as  created  exactly  at  x^  with  a slight 
error.  The  expected  nvimber  of  photons  per  electron  created 
along  the  path  is 


-// 

-'O  -'O 


d rE(s) 


drt' 

(E(s),E'y)  dE/y.  ds 

dEy  • ' 


(1^2) 


where  d^/dEy  is  the  macroscopic  photon  production  cross  section, 

i.e.,  the  probability  per  unit  path  length  and  per  unit  photon 
energy  of  producing  a photon  with  energy  Ey, 

Several  modeling  possibilities  arise.  One  is  to  use  the 
expected  number  of  photons  produced  and  apply  Poisson  statistics 
to  compute  the  actual  nxjraber 


^n  = 


V 


n 


exp(-  2/  ) 


n = 0,  1,  ...  (1|3) 


nl 


The  actual  value  of  n is  obtained  by  random  sampling  this  discrete 
density  function.  Each  photon  is  given  a weight  equal  to  that 
of  the  electron  and  saved  for  later  processing.  At  the  same 
tim.e,  the  electron  energy  can  be  decremented  either  for  each 
photon  so  produced,  or  for  the  average  energy  of  the  produced 
photons . 


If  this  electron  energy  decrementing  procedure  is  used,  the 
value  of  L(E)  used  in  obtaining  the  free-flight  path  must  exclude 


21 


i 


the  radiative  stopping  power.  Otherwise  the  electron  loses  this 
radiative  source  energy  twice. 


An  alternate  procedure  is  to  produce  photons  with  weight 
equal  to  the  expected  production  and  use  the  total  L(E)  with 
no  electron  energy  decrement  during  the  electron  transport. 

Electrons  can  also  excite  the  orbital  electrons  causing 
fluorescent  photons.  The  expected  mamber  of  excitations  along 
the  free-flight  path  is 


Hi  = f Cf^(E 

-'O 


(s))  ds 


(4i+) 


where  ^'^(E)  is  the  cross  section  for  exciting  level  i.  Once 

excited,  the  atom  produces  fluorescence  photons  with  probability 
where  is  the  fractional  fluorescence  yield  given  excita- 
tion. Auger  electrons  and  photoelectrons  are  also  produced. 

All  are  saved  for  later  processing. 

The  electron  interaction  modeling  also  includes  the  production 
of  high  energy  secondary  electrons.  The  expected  number  of 
secondary  electrons  produced  with  energies  greater  than  A is 
computed  as 


nE(s)-A 
dS' 


’ ) dE'  ds 


OtS) 
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where  d(^ /dS'  is  the  macroscopic  cross  section  per  unit  path 
length  and  per  unit  degraded  energy  for  an  electron  to  produce 
secondary  electrons  with  energy  E'  and  E”  where  E = E*  + E” . 

The  cross  section  is  given  by  Holler  for  electron-electron  inter- 
actions. Two  electrons  always  come  out  of  this  reaction.  The 
electron  with  energy  greater  than  E/2  is  considered  the  primary 
and  the  electron  with  energy  less  than  E/2  is  called  the  secondary 
electron. 

Given  the  expected  number  of  high  energy  secondary  electrons, 
there  are  several  ways  to  proceed.  The  procedures  are  similar 
in  logic  to  the  creation  of  bremsstrahlung  photons.  One  model 
uses  Poisson  statistics  to  detemine  the  actual  number  of  created 
secondaries.  Another  model  creates  one  secondary  electron  with 
a relative  weight  equal  to  the  expected  number  of  secondaries. 

In  both  instances,  the  initial  electron  energy  can  be  decremented 
by  the  actual  or  expected  secondary  energy.  If  the  primary 
electron  energy  is  decremented,  then  L(E)  used  to  obtain  the 
free-flight  path  must  exclude  those  energy  losses  which  correspond 
with  the  creation  of  secondary  electrons. 

The  energy  S'  of  any  secondary  electrons  is  determined  by 
sampling  the  spectrum  of  secondary  electron  energies.  This 
spectrvim  is  simply  the  normalized  Holler  cross  section  which 
is  stored  as  a tabulated  cxomulative  distribution  function  on 
a fixed  mesh  of  primary  electron  energies  and  of  secondary 
electron  energies  (expressed  as  a fraction  of  the  primary  energy). 

The  direction  of  secondary  electrons  relative  to  the  primary 
electron  and  the  deflection  of  the  primary  due  to  secondary 
production  are  determined  by  the  energy-angle  relationship  for 
electron-electron  scattering. 
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The  last  model  applied  to  primary  electrons  accounts  for 
energy  loss  fluctuations.  Here  a Gaussian  formula  is  used 


p(E”  - E' ) 


E»  )^  / 2 O'  ^ ] 


(ii6) 


where  O'  is  tlie  mean-square  deviation  of  the  electron  energy 
from  the  mean  for  the  free-flight  path  and  is  computed  as  the 
path  length  integral 

,E(s) 

Cf=  (E(s),E»)  (E-E')^dE«ds  (1^7) 

J J 

0 -'o 

where  d^ /dE’  is  the  differential  energy  loss  cross  section, 
i.e.,  the  Moller  electron-electron  scattering  cross  section. 

The  mean-square  energy  loss  is  stored  as  a tabulated  table  on 
a mesh  of  discrete  energy  points. 

Energy  loss  fluctuations  from  radiative  processes  are  modeled 
explicitly  in  the  treatment  of  bremsstrahlung  production.  The 
secondary  electron  models  which  explicitly  correct  the  primary 
electron  energy  for  the  produced  secondary  electron  must  exclude 
this  energy  loss  range  in  the  computation  of  the  mean-square 
loss,  i.e.,  the  upper  limit  on  the  inner  integration  becomes 
E(s)  - A . 

These  same  models  are  applied  on  each  condensed  history  step 
of  the  electron  until  it  escapes  from  the  geometry  or  until  the 
energy  becomes  so  small  that  the  electron  residual  range  is  much 
less  than  the  dimensions  of  the  material  geometry.  With  small 
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residual  energy,  the  flux  spectrum  as  the  particle  goes  to  zero 
energy  is  determined  from 


) = W E / L(E» ) 


0 ^ E'  * E (46) 


During  the  electron  tracking  process  numerous  scores,  based 
on  the  history  of  the  particle,  can  be  made.  The  obvious  ones 
on  a macroscopic  scale  are  entrance,  transmission,  reflection, 
absorption,  and  energy  deposited  relative  to  the  defined  problem 
geometry. 

From  the  definition  that  collision  density  equals  flux  times 
the  total  cross  section,  the  flux  in  the  slab  can  be  estimated 
by  dividing  the  particle  weights  going  into  collisions  by  the 
total  cross  section  at  their  collision  energy.  For  the  condensed 
history  technique,  the  implied  total  cross  section  is 


(f(E)  = L(E)/f  E 


(49) 


where  f is  the  fractional  energy  loss  per  step.  This  same  flux 
definition  can  be  used  to  calculate  scattered  flux  at  an  arbitrary 
point  since  the  flux  into  collision  times  the  differential 
scattering  cross  section  yields  the  scattered  source  for  all 
directions,  in  particular,  the  direction  towards  a detector  point. 
Finally,  the  particle  weight  coming  out  of  collision  times  the 
collision  probability  can  be  integrated  along  the  path  of  poten- 
tial collisions  to  yield  the  track  length  flux  definition. 
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Photon  Tracking 


When  the  electron  and  any  secondary  electrons  have  been 
completely  tracked,  secondary  photons  are  created  and  tracked. 

Two  modes  of  photon  creation  were  discussed  during  the  electron 
tracking  process.  One  mode  is  bremsstrahlung,  which  is  forward 
peaked  in  angular  distribution  relative  to  the  direction  of  the 
electron.  The  other  mode  is  direct  excitation  fluorescence 
which  is  isotropic. 

Since  every  electron  path  yields  an  expected  number  of  both 
bremsstrahlung  and  direct  excitation  photons,  both  secondary 
production  modes  are  considered.  A decision  is  made  at  random 
at  each  photon  production  site  (terminal  points  of  the  electron 
steps)  to  generate  and  track  either  the  bremsstrahlung  or  the 
fluorescence  phosons.  The  decision  is  based  on  the  total  energy 
of  the  photons  expected  for  the  two  modes  and  the  initial  weight 
of  the  photon  is  adjusted  for  this  decision  process. 

Bremsstrahlung  photons  are  created  with  an  energy  sampled 
from  the  normalized  differential  (per  energy)  production  cross 
section.  The  direction,  position,  and  time  are  set  to  the  values 
of  the  primary  electron  at  the  step  where  the  photon  was  produced. 

This  bremsstrahlung  model  assumes  that  the  photons  are 
produced  in  the  same  direction  as  the  primary  electron.  Alterna- 
tively, the  direction  of  bremsstrahlung  photons  relative  to  the 
electron  is  sampled  from  the  double  differential  bremsstrahlung 
cross  section  normalized  to  form  a probability  density  function. 
This  requires  an  associated  equiprobable  azimuth  selection  and 
a rotation  from  the  electron  direction  reference  frame  to  the 
geometry  reference  frame. 
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Direct  excitation  photons  are  created  with  spectral  weights 
based  on  the  primary  electron  energy,  weight,  and  step  size 


= W d CXiCE) 


i = 1,  2,  ...  (50) 


i.e.,  a photon  for  each  excited  level.  This  spectrum  is  then 
sampled  to  obtain  a single  m onoenergetic  photon.  The  photon 
direction  is  isotropic  and  is  sampled  directly  in  the  geometry 
reference  frame. 

Once  created,  the  photon  is  tracked  through  the  material 
slab  until  it  leaks  out,  degrades  below  energy  cutoff,  or  is 
absorbed. 

The  distance  between  photon  collision  points  is  calculated 
by  one  of  several  alternate  methods.  A simple  particle  analogy 
yields  the  distance  (for  a single  material)  of 


d 


J 

O'(E) 


ln(1 


r) 


(51  ) 


If  the  distance  d exceeds  the  distance  to  the  boundary  of  the 
geometry,  the  photon  has  its  collision  in  some  part  of  the  viniverse 
which  was  not  modeled.  Ey  default,  this  collision  is  not  germane 
to  the  problem.  The  photon  escapes  from  the  geometry  and  its 
history  is  terminated. 

If  the  distance  d implies  a collision  within  the  geometry, 
the  next  step  involves  the  treatment  of  the  interaction  at  the 
selected  collision  point.  Before  examining  the  interaction,  it 
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is  instructive  to  examine  other  methods  for  obtaining  the  distance 
between  photon  collisions. 


In  particular,  the  probability  density  function  for  collision 
points  can  be  normalized  to  yield  a collision  point  within  the 
defined  geometry 

w CT(E)  exp(  -O'(E)  s ) 

p'(s)  = (52) 

1 - exp(  -Cf('S)  s^  ) 

where  s^  is  the  distance  along  the  photon  direction  to  the 

boundary  of  the  material  geometry.  This  probability  density 
forces  collisions  within  the  defined  geometry  but  requires  a 
.particle  weight  adjustment  factor  to  compensate  for  the  biased 
sampling 

p(s) 

W«  = W = W (1  - exp(  - (f{E)  s,  ))  (53) 

P (s) 

i.e.,  the  particle  weight  is  reduced  by  the  nonescape  probability. 

Other  collision  density  probability  fiinctions  are  often  used. 
This  is  a relatively  easy  way  to  modify  the  physics  of  the  sampled 
interactions  to  obtain  those  interactions  deemed  important  to  the 
transport  problem.  Particle  weights  are  always  adjusted  to 
preserve  expected  values  by  multiplying  by  the  ratio  of  the 
true  probability  density  and  the  modified  probability  density. 
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The  photon  interaction  event  is  selected  at  random  by  sampling 
the  discrete  probability  density  formed  by  the  ratio  of  individual 
interaction  cross  sections  to  the  total  cross  section.  For  some 
problems,  this  \inbiased  interaction  selection  is  modified  by  an 
estimated  importance  of  each  interaction  with  a corresponding 
modification  of  the  particle  weight. 

If  Compton  scattering  is  selected,  the  energy  of  the  scattered 
photon  is  obtained  by  a rejection  technique  using 


X = E/S' 


1 1 

g(x)  = — 

ln(1  + 2 T)  ^ 


X 2 1 + 2 T 


h(x)  = X (x  + 1/x  - 1 + ) / 2 


= 1 - (x  - 1)/T 

T = E / m^c^  iSk) 


The  direction  and  energy  of  both  the  photon  and  the  recoil 
electron  are  determined  by  the  mechanics  of  the  Compton  inter- 
action. 


If  incoherent  form  factors  are  available,  the  photon  and 
electron  weights  are  multiplied  by  the  bound  electron  correction 
factor,  which  is  a function  of  element  atomic  number  Z, 


W = W S(q,Z)  / Z 


(55) 


which  is  evaluated  by  linear  interpolation  of  a tabulation 
versus  q (momentum  transferred) 


( 1 - 2^/x  + 1/x^ 


(56) 


If  photons  are  being  tracked  as  polarized  particles,  a 
polarization  vector  e perpendicular  to  u is  assumed  (generated 
at  random  at  the  photon  birth  point).  The  decision  to  follow 
the  polarized  or  unpolarized  photon  out  of  the  interaction  is 
based  on  the  relative  cross  sections.  If  the  perpendicular 
(random  polarization  after  scattering)  component  is  picked,  an 
azimuthal  angle  is  determined  from  the  equiprobable  azimuthal 
distribution.  If  the  parallel  mode  is  selected,  the  azimuthal 
angle  is  also  sampled  from  the  equiprobable  distribution  and 
the  photon  weight  is  adjusted  to  correct  for  this  biased  sampling. 
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If  coherent  scattering  is  the  selected  photon  interaction, 
the  cosine  of  the  scattered  angle  is  selected  from  the  angular 
probability  density  fxmction  using  a rejection  technique 


p(M^  = 3 (1  + / 8 


g(/i)  =1/2 


h(/i)  = (1  + / 2 


(57) 


The  weight  of  the  photon  is  then  reduced  by  the  form  factor 


j 

\ 


\ 


W = W F^(q,Z)  / 


(58) 


which  is  tabulated  versus  momentum,  q,  transferred 

q = —V  (2  (1  - /I)  ) ^ (59) 

m c 
o 

If  the  photon  was  polarized,  its  azimuthal  angle  relative 
to  the  polarization  vector  e is  left  unchanged.  If  the  photon 
is  unpolarized,  a random  azimuthal  angle  is  generated.  The 
resulting  u'  direction  relative  to  u.  is  then  rotated  to  the 
reference  coordinate  frame. 


1 
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If  pair  production  is  the  selected  photon  interaction,  two 
annihilation  photons  are  generated  at  the  interaction  site  by 
increasing  the  photon  weight  by  two,  setting  the  energy  to 
and  generating  an  isotropic  direction. 

This  pair  procedure  is  modified  if  pulse  height  distributions 
are  being  obtained,  since  it  is  necessary  to  track  each  of  the 
annihilation  photons  separately. 

The  positron/electron  pair  are  created  by  sampling  an  equi- 
probable  energy  distribution  on  (0,E  - 2m^c  ).  Both  members  of 
the  pair  are  tracked  as  electrons.  Their  direction  is  selected 
from  an  exponential  distribution  which  preserves  the  mean 
production  direction. 

If  the  photoelectric  event  is  the  selected  interaction,  the 
subshell  of  the  interaction  is  determined  proportional  to  the 
partial  cross  section  ^^(E).  The  photon  weight  is  multiplied 
by  the  fluorescence  yield  and  the  photon  energy  is  set  to  the 
fluorescence  energy.  The  corresponding  direction  is  isotropic 
and  is  selected  at  random  in  the  reference  coordinate  frame. 


The  Auger  electrons  associated  with  the  photoelectric  event 
are  given  isotropic  directions.  The  photo  electron  direction 
relative  to  the  photon  direction  is  sampled  using  the  rejection 
technique 


n"  - ' 


2 


(7?  - 


h(/X)  = (7J^  - 1)  (1  . 

^ = m^c2  / (E  (s  + 2 m^c2  ) ) (^qj 

Ahe  weight  of  the  photo  electron  is  multiplied  by  the  factor 

p /p''(/z)  where  p (M)  i3  the  Sauter  or  Fischer  angular 

distribution  and 


) = g(//)  h(//)  / G 


(1  - jl^) 
(7l  - 


(61) 


The  photon  tracking  continues  by  collision  using  similar 
logic  on  each  collision  until  the  photon  escapes,  its  energy 
goes  below  a cutoff,  or  its  weight  becomes  small. 

During  the  photon  tracking,  scoring  comparable  to  that  for 
the  primary  electrons  is  performed. 
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SECTION  III 
ANALYSIS  PROCEDURES 

The  MASTER  program  file  contains  several  particle  tracking 
procedures.  Most  of  these  procedures  utilize  part  or  all  of 
the  geometry,  source,  and  detector  data  base  implied  in  the 
previous  section  and  discussed  in  more  detail  in  the  next  section. 

BETA-II  Procedure  (refs.  5 and  6) 

BETA  performs  Monte  Carlo  electron-bremsstrahliong  or  photon 
electron  transport  by  the  methods  outlined  in  Section  II,  pre- 
dominantly analog  Monte  Carlo. 

BETA  uses  the  complete  data  base.  Particle  tracking  is 
controlled  by  subprogram  SOBER  and  includes  an  optional  ability 
of  treating  the  effect  of  fixed  electric  and  magnetic  fields. 

The  fixed  electric  field  option  includes  the  nonlinear  buildup 
of  electric  fields  assuming  frozen  charge.  This  charge  buildup 
is  limited  to  two-dimensional  cylindrical  geometries.  The 
solution  of  Poisson’s  equation  for  a fixed  charge  density  is 
performed  by  subprogram  VOLTS  in  the  output  module  of  the  program. 

For  problems  not  requiring  field  effects,  the  EEAI^A?  procedure 
is  recommended  for  its  superior  particle  processing  speed. 
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BEAI'IAP  Procedure 

BEAMAP  is  a faster  version  of  BETA  where  no  field  effects 
are  allowed  (except  for  an  electron  reflection  boimdary)  and  all 
of  the  appropriate  particle  physics  are  collected  into  a single 


TEMPER  Procedure  (refs.  7 and  6) 
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TEMPER  calculates  time-dependent  charge  and  current  densities, 
scalar  and  vector  potentials,  and  electric  and  magnetic  fields 
in  three-dimensional  cavity  geometries.  A problem  is  initiated 
from  a specified  electron  emission  distribution.  The  time 
evolution  of  a problem  is  obtained  by  a series  of  time  steps  in 
which  potentials  are  assumed  constant.  During  each  time  step, 
typical  electrons  are  sampled  from  the  electron  emission  distri- 
butions and  added  to  the  electron  population  which  existed  at 
the  end  of  the  previous  time  step.  Electrons  are  moved  using 
relativistic  kinematics  which  include  Lorentz  force,  slowing 
down  from  gas  in  the  cavity,  and  production  of  primary  ioniza- 
tion electrons.  At  the  end  of  a time  step,  a snapshot  is  taken 
of  the  time-averaged  charge  and  current  densities  during  that 
step.  At  the  end  of  each  time  step,  the  scalar  and  vector 
potentials  are  calculated  by  performing  voliime  integrations  over 
the  charge  and  current  densities.  Retarded  time  effects  are 
explicitly  calculated. 

I 

DIODE  Procedure 

This  analysis  procedure  calculates  quasi-static  nonlinear 
electron  transport  between  parallel  plates.  Electron  emission 
is  allowed  from  both  plates.  The  energy  dependence  of  emission 
is  represented  within  a fixed  group  structure  where  each  group 
can  be  subdivided  for  better  representation  of  highly  charge 
limited  problems. 

The  angular  dependence  of  emission  must  be  an  integral 
power  of  the  cosine  of  the  emission  angle  relative  to  the  cavity 
surface  nomal.  This  angular  dependence  is  represented  by  a 
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set  of  discrete  intervals  which  preserve  the  total  emission  and 
the  average  emission  angle. 

Time  dependence  of  the  emission  is  either  triangular  or 
Poisson,  The  time  dependent  behavior  of  the  cavity  is  determined 
by  solving  a series  of  quasi-static  problems  at  successive  points 
in  time. 


The  cavity  length  is  subdivided  into  intervals  for  representa- 
tion of  the  spatial  behavior  of  the  electric  field  and  of  the 
electron  transport  between  the  plates.  At  each  time  step  typical 
electrons  are  started  from  each  emission  surface  and  in  every 
energy  interval  and  direction  interval.  The  contribution  from 
these  typical  electrons  to  the  charge  density,  current  density, 
and  cavity  gas  ionization  are  computed.  The  charge  density 
and  the  net  current  density,  in  conjunction  with  a fixed  impe- 
dence  between  the  plates,  are  used  to  obtain  the  potential  and 
electric  field, 

EMP3D  Procedure 


This  analysis  procedure  calculates  nonlinear  electron  trans- 
port in  regular  three  dimensional  geometries.  A finite  difference 
particle-in-cell  approach  is  used. 

Particle  populations  are  represented  by  a fixed  mesh  in 
space,  energy,  and  direction.  Emission  distributions  have  the 
sane  characteristics  and  can  be  defined  at  any  mesh  point. 
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This  analysis  procedure  calculates  neutron/photon  transport 
using  Monte  Carlo  methods.  Extensive  importance  sampling  is 
used. 


SOBER  uses  the  complete  data  base  to  obtain  volume  and 
surface  averaged  fluxes.  Fluxes  are  also  obtained  at  multiple 
point  detectors.  In  general,  this  procedure  is  not  used  for 
point  detectors  in  source  and  scattering  volumes  where  explicit 
representation  of  point-to-point  singularities  are  required. 

The  PASTER  analysis  procedure  is  used  for  these  point  detector 
problems. 

PASTER  Procedure  (refs.  1-3) 

FASTER  is  used  to  obtain  f lux-at-a-point  when  the  point  is 
in  a source  or  scattering  volume.  Monte  Carlo  methods  are  used 
to  obtain  both  neutron  or  photon  fluxes. 


This  procedure  uses  importance  sampling  which  includes  the 
point  source  and  point  detector,  singularities . Therefore,  finite 
variance  fluxes  are  always  obtained. 


The  complete  data  base  is  used  by  FASTER  including  options  for 
transport  in  helical  ducts  (ref.  12)  and  shield  weight  optimiza- 
tion (ref.  13). 


37 


PULSE,  TRICK,  TREAT  Procedures  (refs.  and  15) 


PULSE  performs  parametric  electron  transport  calculations  of 
charge  and  energy  deposition,  reflection,  transmission,  and  pulse 
height  distributions.  The  procedure  combines  Monte  Carlo  and 
numerical  integration  procedures.  An  electron  track  is  initiated 
at  the  maximum  energy  of  interest  and  followed  through  an  infinite 
medium  by  Monte  Carlo  methods.  The  track  at  various  degraded 
energies  is  used  to  characterize  electron  sources  at  these  lower 
energies.  The  track,  or  any  residual  part  used  for  lower  energies, 
is  then  rotated  and  translated  nximerically  to  yield  information 
for  various  initial  directions  and  material  thicknesses.  The 
procedure  is  repeated  for  multiple  tracks  to  converge  the  Monte 
Carlo  aspect. 

TRICK  and  TREAT  calculate  pulse  height  distributions  for 
small  volume  detectors  inside  complicated  geometries  by  mmieri- 
cal  integration  using  either  approximate  transport  kernels  (TRICK) 
or  interpolations  of  tabulated  Monte  Carlo  transport  data  (TREAT). 
Both  prima.'y  and  secondary  radiations  contribute  to  i;he  pulse 
heights.  Major  allowed  primary  sources  are  neutrons  producing 
activation,  photons,  and  electrons. 

AKALOG  Procedure 

This  analysis  procedure  uses  analog  Monte  Carlo  techniques 
to  simulate  photon  interactions  within  a complex  geometry 
pulse  height  detector  system.  The  detector  is  assumed  to  reside 
within  a laboratory  mockup.  Importance  sampling  is  used  to 
increase  the  number  of  particles  which  enter  the  detector  both 
directly  or  by  scattering  from  the  facility  structure. 


SHIELD  and  SIGMA  Procedures  (refs.  16  and  1?) 


SHIELD  numerically  integrates  the  electron-bremsstrahlung 
transport  equation  for  one  dimensional  geometries  to  obtain 
dose  and  other  fl;ix  responses  versus  shield  thickness.  Cross 
section  data  and  response  functions  are  obtained  from  the  data 
base . 

SIGMA  calculates  space  radiation  dose  at  points  in  compli- 
cated geometries  by  a solid  angle  integration  about  the  point 
using  interpolated  one-dimensional  dose  transmission  kernels, 
e.g.  SlilELD  results.  SIGMA  uses  the  geometry  and  detector 
portions  of  the  data  base. 

ADJOINT  Procedure  (refs.  l8  and  19) 

This  Monte  Carlo  procedure  performs  electron  and  photon 
transport  calculations  by  tracking  particles  backward  from  the 
detector  to  the  source.  It  is  used  for  point  or  small  detector 
calculations  in  extended  source  fields,  e.g.,  space  radiation 
sources.  Capabilities  include  dose  calculations,  pulse  height 
distributions,  and  parametric  dose  attenuation  data  for  electron 
sources  and  secondary  electron  flux/dose  for  primary  photon 
sources.  The  complete  data  base  is  used. 


SECTION  IV 
DATA  BASE 

The  data  base  is  shared  between  subprograms  through  named 
common  blocks  and  blank  common.  The  layout  of  these  areas  is 

a.  Blank  common  contains  data  arrays.  Each  data  array 
has  variable  dimensions  and  arrays  are  packed  into 

ascending  locations  according  to  the  current  array  i| 

dimensions . 

b.  Named  common  blocks  contain  fixed  location  data.  In 
particular,  the  addresses  and  dimensions  of  data  arrays 
are  stored  in  named  common  blocks. 

These  areas  are  initialized  before  processing  user  data  by 
subprograms  EOQTIT,  which  initializes  executive  data  arrays, 

DATBAS,  which  allocates  arrays  used  for  problem  description, 
and  DATHUN,  which  allocates  arrays  used  during  the  problem 
solution,  II 

Input,  preparation,  storage,  and  use  of  the  data  base  is 
divided  into  many  procedures.  This  division  simplifies  the 
description  and  modification  of  transport  problems.  Each  pro- 
cedure is  controlled  by  user  directives.  These  procedures  are 
described  under  the  general  categories 

a.  Executive  procedures 

b.  Geometry  | 

c.  Source  ; 

d.  Detector 

e.  Physics 
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Executive  Procedures 


All  procedures  are  controlled  by  user  supplied  data.  Each 
procedure  involves  a block  nf  consecutive  data  cards.  The  first 
card  of  the  block  is  the  "header  card”  and  contains  the  name  of 
the  procedure,  where  the  first  three  letters  are  unique,  e.g., 

STOp.  The  second  card  of  the  block  is  the  "lirait/option  card" 
and  specifies,  as  integer  data,  array  limits  and  options  for 
the  designated  procedure. 

'wTien  a procedure  is  recognized,  the  limits  and  options  are 
stored  in  appropriate  locations  of  the  named  common  blocks. 

Array  storage  is  reallocated  (subprogram  lARHAY)  if  any  array 
dimensions  have  changed.  If  the  procedure  requires  local  data 
arrays,  these  arrays  are  allocated  to  blank  common  at  the  end 
of  the  data  base  arrays  using  subprogram  IIIDEX.  This  local 
data  storage  area  is  recovered  after  the  procedure  is  completed. 

Input  data  required  by  a procedure  is  processed  by  subprogram 
READER.  During  the  input  process,  the  following  options  are 
available 

a.  Fixed  or  variable  format  data  with  variable  format 
data  separated  by  blanks  or  commas 

b.  Repeat  and  interpolation  of  data  elements.  Interpola- 
tions include  exponential,  logarithmic,  linear  and 
power  laws  2 through  9 

c.  Conversion  of  specified  data  elements,  e.g.,  degrees  to 

radians,  exponentiation,  and  other  function  evaluations  - , 
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A problem  is  terminated  by  the  STOp  procedure.  An  optional  for- 
matted dump  of  the  data  base  is  obtained  prior  to  the  programmed 
stop.  This  data  base  dump  includes  each  named  common  block  and 
each  variable  dimension  array. 


A dump  of  the  entire  data  base,  or  of  specified  named  commons 
and  arrays,  can  also  be  obtained  at  the  entrance  and  exit  points 
of  major  subprograms.  The  necessary  control  information  is 
supplied  as  data  using  the  DUMp  procedure. 


Arguments  passed  to  the  data  card  processor  indicate  whether 
the  input  data  has  specific  units  and  order.  This  information 
permits  data  conversion  to  standard  units  during  the  input  pro- 
cess. It  also  allows  checking  arrays  for  order  consistency. 

Other  errors  detected  during  the  input  process  include  data  out 
of  range,  indices  out  of  range,  and  missing  data. 

Obvious  executive  procedures  are  the  call,  by  name,  of  the 
analysis  procedures,  e.g.,  BEAmap.  For  analysis  procedures, 
the  limit/option  card  contains  the  number  of  histories  and  other 
integers  which  control  the  particle  tracking.  When  an  analysis 
procedure  is  recognized,  the  following  operations  are  performed 

a.  Missing  data  is  noted  and  where  possible  supplied  by 
the  program 

b.  The  dimensions  of  data  arrays  used  for  particle  tracking 
are  set  to  necessary  values  for  problem  solution 

c.  If  the  problem  setup  has  no  detected  errors,  the  analysis 
is  performed.  If  errors  were  detected,  the  analysis  step 
is  aborted. 


Long  rvinning  problems  can  be  solved  using  multiple  runs. 

A restart  file  is  specified  using  the  TAPes  procedure.  The 
information  written  on  the  restart  file  is  retrieved  in  a subse- 
quent computer  run  using  the  COKtinue  procedure.  The  analysis 
is  continued  using  the  history  counter,  accumulated  scoring, 
etc.  contained  on  the  restart  file. 

Trivial  executive  procedures  include 

a.  LABel,  for  entering  printout  headings 

b.  HSXt,  to  increment  the  problem  identification  number 

c.  ZSRout,  to  re-initialize  between  independent  problems 

d.  RAl'idom,  to  override  the  default  random  number 

Some  executive  procedures  alter,  in  some  manner,  the  pro- 
cessing of  data  by  the  problem  related  processors 

a.  PRInt,  suppresses  the  printout  of  extensive  input  cards 

b.  UNIts,  indicates  that  some  inputs  have  nonstandard 
physical  \inits,  e.g.,  energies  in  keV,  geometry  dimen- 
sions in  feet,  etc. 

c.  ADDress,  supplies  index  modifiers  for  problem  data 
input  by  index,  e.g.  geometry,  sources,  and  detectors. 
This  procedure  simplifies  the  collation  of  two  indepen- 
dent problems  into  a single  problem. 

Finally,  data  can  be  directed  to  the  data  base  using 

a.  BLOck,  specifies  the  value  of  elements  of  named  common 
blocks 

b.  ARRay,  accepts  variable  dimension  arrays  by  name 
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Geometry 

The  transport  problem  geometry  is  described  by  disjoint 
voliames  called  ’’regions"  where  the  boxindaries  of  regions  are 
called  "surfaces".  Each  region  has  constant  material  properties 
denoted  by  a composition  or  material  identifier  and  a density 
relative  to  the  reference  density  of  that  material. 

The  problem  geometry  can  be  compact,  i.e.,  containing  no 
undefined  voliomes,  or  contain  undefined  volumes.  All  undefined 
volvimes  are  assumed  to  be  void.  In  particular,  the  space  sur- 
rounding a compact  geometry  is  assumed  to  be  void.  A compact 
geometry  implies  that  the  surface  forming  the  boundary  betvjeen 
adjacent  regions  is  unique. 

Geometry  related  calculations  assume  that  all  surfaces  are 
described  by  the  quadric  equation 


U(x)  = Sq  + a^  X 1 + a^  z 

2 2 2 
+ a^  X + a^  y + z 

+ a^  xy  + a^  yz  + a^  zx  (62) 

where  a^^  , i = 0,  1,  ...,  9,  are  constants.  U(x)  is  zero  for 

points  on  the  surface,  U(x)  is  less  than  zero  for  points  "inside" 

the  surface,  and  U(x)  is  greater  than  zero  for  points  "outside" 

the  surface.  The  inside/outside  sense  can  be  reversed  by  multi- 
plying the  surface  equation  by  -1  . 


Description  of  a geometry  using  these  quadric  surfaces 
is  tedious.  Therefore,  a number  of  processors  are  used  to 
simplify  the  description  task.  These  processors,  as  a group, 
are  contained  in  subprogram  GEOMIIj. 

Geometry  related  calculations  include 

a.  Determining  the  region  occupied  by  a point,  performed 
by  subprogram  LOCATE 

b.  Ray  tracing  from  a point  in  a fixed  direction  to  deter- 
mine the  order  and  thickness  of  materials  encountered 
along  a particle  path,  performed  by  PATH 

c.  Calculating  the  surface  normal  at  specific  boundary 
crossings,  performed  by  subprogram  NORl-IAL 

The  numerical  operations  performed  during  these  calculations  are 
detailed  in  references  1,  S,  and  20. 

Some  transport  problems  have  simple  geometries.  It  is 
sxifficient  for  these  problems  to  describe  the  material  geometry 
using  one  of  the  following  processors 

a.  SPHere,  for  concentric  spherical  zones.  For  the  parti- 
cular problem  of  transport  in  an  approximation  of  a 
spherical  exponential  atmosphere,  the  AIR  processor 
supplies  the  density  versus  altitude  variation. 

b.  CYLinder,  for  cylindrical  geometries  with  azimuthal 
symmetry 

c.  SIHple,  for  regular  meshing  of  the  three  orthogonal 
coordinates  of  rectangular,  cylindrical,  or  spherical 
geometries 

These  procedures  assume  a single  material  in  all  regions.  This 
assumption  is  countermanded  using  tne  ARRay  processor  for  the 
specific  array  MIL  which  contains  the  material  identifier. 
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More  complicated  geometries  often  use  direct  description 
of  surfaces  and  regions.  These  descriptions  can  be  used  to 
augment  a subset  of  the  geometry  described  by  the  more  simple 
geometries,  or  these  descriptions  can  be  used  to  specify  the 
entire  geometry. 

Direct  surface  description  is  accomplished  with  the  SURface 
processor.  This  processor  accepts  the  coefficients  for  the 
general  quadric  surface.  Alternatively,  simple  surfaces  are 
recognized  in  more  simple  forms  and  expanded  by  the  processor 
to  the  quadric  form.  These  simple  forms  include 

a.  Planes  perpendicular  to  one  of  the  coordinate  axes 

b.  Planes  parallel  to  one  of  the  coordinate  axes 

c.  Cones  parallel  to  a coordinate  axis 

d.  Cylinders  and  elliptic  cylinders  parallel  to  an  axis 

e.  Spheres  and  ellipsoids 

Some  regions  have  helicoid  or  torroidal  boundaries.  The  cor- 
responding surface  equations  are  specified  with  the  HELix  pro- 
cessor. These  surface  equations  utilize  a special  intersection 
calculation  during  ray  tracing  (ref.  21). 

Direct  region  description  is  performed  via  the  REGion 
processor.  The  description  includes 

a.  Material  identifier  and  relative  density 

b.  Bovindary  surface  list 

c.  Coordinates  of  any  point  in  the  region 

The  point-in-region  coordinates  are  used  to  determine  the  in- 
side/outside sense  of  the  region  relative  to  the  boundary  sur- 
faces . 
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A SURface  and  REGion  description  of  problem  geometries  can 
also  get  tedious.  This  has  led  to  the  following  processors, 
which  have  a limited  applicability 

a.  BOArds,  for  describing  parallel  material  slabs  with 
different  compositions  and  common  transverse  boundaries 

b.  BOXes,  for  describing  multiple  layer  boxes,  cylinders, 
and  spheres  embedded  in  regions 

c.  BAYS,  for  describing  the  exterior  skin  and  interfaces 
between  multiple  bay  spacecraft 

d.  Divide,  for  subdividing  regions  into  smaller  components 

A more  general  procedure  for  simple  description  of  regular 
geometry  regions  is  DSSign.  This  procedure  combines  the  surface 
and  region  description  for  the  following  shapes 

a.  plates 

b.  elliptical  cylinders 

c.  spheres  and  ellipsoids 

d.  cylindrical  annuli 

e.  truncated  cones 


This  procedure  shortens  the  geometry  description  but  requires 
extreme  care  in  data  preparation  so  that  surfaces  forming  a 
common  boundary  between  two  regions  are  recognized,  i.e.,  dim- 
ensions supplied  with  exactly  the  same  significant  digits.  The 
procedure  assists  in  this  recognition  by  checking  for  voliwietric 
overlap  of  the  regions. 

Some  geometries  involve  regular  surfaces  not  aligned  with 
the  coordinate  axes.  The  description  of  these  surfaces  is 
accomodated  directly  with  the  SURface  processor.  The  descrip- 
tion is  simplified,  however,  by  using  the  ROTate  processor. 


The  ROTate  processor  operates  on  specified  surfaces  and 
regions.  Operations  include  rotations  or  translations,  any 
number  and  in  any  order.  The  same  processor  can  be  applied  to 
the  geometry  of  fixed  sources  and  detectors. 

The  description  of  some  problem  geometries  is  simplified 
by  invoking  the  SYl'Imetry  processor  to  delineate  the  siirfaces 
of  symmetry.  However,  this  procedure  cannot  be  used  for  analysis 
procediires  which  perform  f lux-at-a-point  calculations. 
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Verification  of  a geometry  description  includes  several 
self-consistency  checks  performed  by  subprogram  REGION.  The 
most  useful  checking  is  provided  by  exercising  the  ray  tracing 
operations  during  the  generation  of  geometry  pictures  using 
the  Picture  processor  for  printout  display  and  the  CAMera 
processor  for  cathode  ray  tube  (CRT)  display.  Both  processors 
require  specification  of  the  geometry  frame  to  be  viewed.  Both 
will  also  project  the  trajectories  of  particles  generated  during 
short  runs  of  the  Monte  Carlo  analysis  procedures. 


Source 


The  fixed  radiation  source  can  be  supplied  with  a generality 
coimnensurate  with  that  of  the  problem  geometry.  Single  or 
multiple  sources  can  be  defined  using  the  SOUrce  processor.  The 
description  of  each  source  includes  the  spatial,  angular,  and 
energy  distribution  of  emitted  particles.  Time  dependence  is 
not  required  for  msiny  problems  and  is  discussed  as  a separable 
topic. 

The  geometric  characteristics  of  each  source  are  specified 
by  the  following  information 

a,  source  geometry,  rectangular,  cylindrical,  or  spherical 

b,  separable  tabulated  distribution  functions  for  each  of 
the  three  spatial  variables 

c,  a translation  vector 

Irregular  shaped  sources  can  be  specified  by  rejecting  regular 
geometry  source  points  which  are  outside  specific  geometric 
regions. 

The  tabulated  spatial  distributions  can  be  a discrete 
point  (a  delta  fianction)  or  multiple  points.  This  flexibility 
yields  point,  line,  surface,  or  volume  distributed  sources. 

The  distribution  functions  are  accepted  as  \mnormalized  relative 
density  fiinctions. 


The  angular  character  of  a source  is  specified  by  a tabulated 
distribution  in  each  of  the  two  angular  variables.  The  tabulated 
distribution  is  an  unnormalized  relative  density  function  with 
a single  point  (delta  function)  or  multiple  points.  This  yields 
mcr.odirection  or  angularly  distributed,  including  isotropic. 
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sources,  depending  upon  the  number  of  points  in  each  distribution. 
The  interrelationships  of  the  spatial  and  angular  variables  is 
discussed  in  the  Users  Guide,  The  nvimerical  techniques  for 
normalizing  and  interpolating  the  tabulated  distributions  are 
given  in  references  1 and  5. 

Many  options  are  provided  for  describing  the  particle  energy 
spectrum  including 

a.  Tabulated  differential  nvunber  spectrum 

b.  Tabulated  differential  intensity  spectrum 

c.  Tabulated  interval  number  or  intensity  spectrum 

d.  Tabulated  integral  nvimber  or  intensity  spectrum 

e.  Analytic  functions  for  fission  neutrons,  black  bodies, 
fission  gammas,  Gaussian,  Poisson,  and  self-rectified 
electron  beam 

The  fixed  source  is  normalized  to  a specified  total  yield 
in  particles  or  energy.  Alternatively,  the  input  source  can 
be  arbitrarily  scaled  or  forced  to  match  a specified  density 
at  a particular  point. 

Time  dependent  problems  with  nonlinear  aspects  require  that 
a time  profile  be  supplied  using  the  PROfile  processor.  For 
linear  problems,  this  time  profile  is  used  in  a convolution  of 
time  de^^endent  response,  i.e.,  in  linear  problems  all  particles 
are  born  at  t = 0.  Extremely  short  time  analyses  of  nonlinear 
problems  require  the  WAVefront  processor  to  specify  the  time  of 
arrival  of  the  radiation  wavefront. 


In  the  absence  of  importance  sampling,  particles  are  sampled 
from  the  fixed  source  in  an  analog  fashion.  The  analog  sampling 
can  be  modified  in  several  ways.  Stratification  by  source,  spatial 
interval,  and  angular  interval  can  be  selected  through  the  option 
card  for  the  analysis  procedure.  The  relative  frequency  of 
sampling  particles  from  each  source  can  be  modified  through  the 
RELative  processor.  The  frequency  distribution  for  each  of  the 
spatial  and  angular  source  variables  can  be  modified  using  the 
RATio  processor.  Finally,  for  volimie  distributed  sources  and 
point  detectors,  the  1/r  and  exponential  fall  off  of  source 
importance  with  separation  distance  can  be  included  by  the  PSSudo 
source  processor. 


Detector 

The  description  of  detectors  has  a similarity  to  the  descrip- 
tion of  sources.  However,  detectors  are  assumed  to  have  the 
same  sensitivity  for  all  positions  and  directions  in  their  sensi- 
tive volume. 

Macroscopic  information  on  particle  transport  is  obtained 
with  the  DSPosition  and  LSAkage  processors.  These  processors 
yield  information  on  the  total  number  and  energy  statistics  and 
no  information  on  variation  with  particle  energy.  For  problems 
with  extensive  importance  sampling,  e.g.,  f lux-at-a-point  simula- 
tion, the  resulting  output  may  have  little,  if  any,  usefulness. 

Similar,  but  more  definitive  detector  information  is  obtained 
from  the  CHAnnel  processor.  Here  the  deposition  in  specified 
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geometric  regions  is  partitioned  into  multichannel  output  of 
deposition  or  energy  loss  per  particle.  This  pulse  height 
information  can  be  convoluted  with  a specified  smearing  distri 
bution  using  the  NOIse  processor. 

More  detailed  transport  results  are  obtained  for  detectors 
described  using  the  DSTector  processor.  These  detectors  can 
consist  of  points,  geometry  regions,  and  boundaries  of  geometry 
regions.  Results  for  several  detectors  can  be  added  using  the 
SAi^Ie  detector  processor. 

Angular  information  for  these  specific  detectors  can  be 
obtained  using  the  A2JGular  processor  which  yields  Legendre 
moments  of  the  angular  flux,  in  particular  the  forward  and 
backward  current.  These  moments  are  reconstructed  to  yield 
the  azimuthally  averaged  angular  flux  distribution.  Angular 
flvLX  information  can  also  be  obtained  from  the  SOLid  angle 
processor  for  specified  solid  angle  intervals. 

Energy  information  for  the  specific  detectors  is  given  in 
the  energy  mesh  used  to  tabulate  cross  sections.  The  FLUx 
processor  can  be  used  to  establish  a different  output  structure. 
The  RESponse  processor  provides  for  weighting  the  energy  depen- 
dent fluxes  with  either  user  specified  response  functions  or 
response  functions  based  on  cross  section  information,  e.g., 
energy  absorbtion  coefficients. 

Time  dependence  for  specific  detectors  is  made  relative 
to  the  minimum  time  of  arrival  using  the  TRAnslate  processor. 

The  time  dependence  is  obtained  by  time  interval  using  the  Tllle 
processor  and  by  an  analytic  moments,  with  reconstruction, 
method  using  the  i-IOMents  processor.  Both  of  these  time  dependent 
outputs  is  convoluted  with  any  source  time  FHCfiles. 
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Sensitivity  of  the  detector  output  can  be  obtained  for 
source  GROups,  EIRth  sources.  Boundaries  crossed,  ORDer  of 
scattering,  and  SCAttering  regions  by  using  the  indicated  pro- 
cessor. 

More  detailed  sensitivities  can  be  obtained  using  the 
CORrelated  px’ocessor.  The  correlation  is  on  any  combination 
of  source  and  i ->.gion  material.  A correlation  for  electron 
transport  problems  can  be  made  only  for  the  void/nonvoid 
configurations. 


The  FASTER  f lux-at-a-point  analysis  procedxire  can  determine 
the  variation  in  flux  levels  with  respect  to  changes  in  region 
dimensions  by  using  the  NORmal  processor  to  request  derivative 
information,  and  the  THIckness  processor  to  specify  incremental 
material  thicknesses.  The  derivative  information  can  also  be 
used  to  determine  a weight  optimized  shield  configuration  using 
the  MlUimvmi  processor  (ref.  13)« 


Importance  sampling  relative  to  detector  response  can  be 
requested  using  the  PREferred  processor  to  specify  a region 
where  particle  interaction  densities  should  be  increased.  Con- 
versely, the  SKOrt  circuit  processor  can  be  used  to  circumvent 
many  of  the  default  importance  sampling  models.  The  OPTimijm 
processor  can  be  used  to  request  estimates  of  better  importance 
sampling  parameters  based  on  variance  partial  derivatives. 


Many  of  the  detector  outputs  can  be  print  plot  displayed 
through  the  QUIck  plot  processor.  The  PLOt  processor  vfill 
display  the  differential  flux  spectrum  if  a drum  plotter  is 
available . 
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Particle  Physics 


Particle  physics  data  consists  predominantly  of  processed 
cross  sections  obtained  from  detailed  data  libraries. 

The  MATerials  processor  supplies  partial  densities  of  elements 
or  isotopes  in  materials  or  compounds.  The  Mixture  processor 
can  be  used  for  neutral  particle  problems  to  specify  materials 
in  which  only  the  hydrogen  content  varies. 

Nonlinear  problems  require  the  DISlectric  processor  to 
specify  macroscopic  electromagnetic  properties.  The  CUHrent 
processor  is  used  to  define  initial  fixed  fields  for  the  BETA 
analysis  procedure. 

Importance  sampling  of  the  particle  physics  is  optional 
and  inc-ludes  BIAsing  to  set  models  for  selecting  particle  collision 
points  and  scattered  directions,  CAPture  to  force  an  analog  simu- 
lation of  neutron  slowing  down,  SKIn  to  artificially  increase 
photon  interactions  near  an  electron  emission  surface,  and 
SPAtial  to  increase  or  decrease  collision  densities  in  material 
regions . 

The  PKOton  processor  prepares  photon  cross-section  data. 

Card  input  is  accepted.  However,  the  usual  source  of  data  is 
Sl'iDF/3  files  23  and  27  (ref.  22). 

The  ELEctron  processor  prepares  electron  cross-section  data. 
Default  models  are  those  of  the  BETA  and  BETA-II  programs 
described  in  refsrences  t and  5.  Electron  brem.s strahlung  pro- 
duction models  from  references  23  and  2i+  are  used.  Alternatively, 
numerical  data  from  reference  25  oan  be  used  for  energies  below 
0.5  MeV.  Direct  excitation  cross-sections  are  calculated  by 
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the  formalism  of  reference  26  but  with  explicit  evaluation  of 
subshell  cross-sections  from  ElIDP/B  file  23. 

The  2L2ctron  processor  will  also  accept  electron/bremsstrah- 
lung  cross-section  data  from  DATAPAC  (refs,  27  and  2:8).  This 
data  is  superior  to  the  default  data  except  for  the  photon  trans- 
port and  direct  excitation  cross-sections.  These  portions  of 
the  default  data  are  used  in  preference  to  the  corresponding 
portions  of  the  DATAPAC  library. 

Heavy  charged  particle  cross-sections  are  also  developed 
by  the  ELEctron  processor  by  selecting  a corresponding  option 
on  the  limit/option  card, 

Keutron  and  secondary  gamma  ray  production  data  are  obtained 
from  El^DF/B  data  files  (ref.  29)  using  the  NEUtron  processor. 
Alternatively,  neutrons  are  simulated  in  the  multigroup  approxi- 
mation using  the  MULtigroup  processor  and  standard  multigroup 
libraries  (ref.  30),  The  multigroup  libraries  usually  include 
cou.-led  photon  transport  data.  The  secondary  photon  production 
data  are  stripped  from  these  libraries  using  the  SECondary  pro- 
cedure, and  the  photon  transport  data  are  discarded  since  the 
point  value  EllDF/3  data  are  preferred. 
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SECTION  V 


RECOMMENDATION  S 

The  generality  of  the  MASTER  program  file  makes  it  difficult 
for  novice  users  to  select  options  and  provide  data  for  a problem. 
The  following  program  modifications  are  recommended  to  eliminate 
the  difficulty 

a.  Provide  complete  problem  descriptions  as  defaults  for 
each  analysis  procedure 

b.  Allow  user  override  of  the  default  data  by  keywords 

For  example,  if  a user  wants  to  solve  the  default  electron- 
bremsstrahlTong  transport  problem,  the  input  data  for  the  default 
problem  could  consist  of 

,EEA1'IAP,ELECTR0N,PH0T0K/  Default  Problem  Execution 

where  the  default  problem  parameters  describe  a tungsten  target, 
plane-parallel  monoenergetic  source  of  1 MeV  electrons  normal 
to  the  target,  and  a detector  measuring  backward  and  forward 
photon  emission. 

To  change  default  parameters,  keyword  parameters  could  include 

, BEAl-IAP , ELECTRON , PHOTON 
, TARGET=GOLD 

, SOURCE. DIRECT ION =30DEGREES 
, SOURCE . ENERGY=1  OOKEV/ 


Similar  default  problems  and  keyword  parameters  would  be  available 
for  each  analysis  processor. 
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For  more  complicated  problem  descriptions  it  is  recommended 

that 

a.  Materials  contained  in  geometric  voliames  be  specified 
by  name  and  density 

b.  Material  compositions  and  standard  densities  be  obtained 
from  a library  that  contains  all  elements  and  the  more 
common  mixtures,  alloys,  and  compounds 

c.  All  macroscopic  particle  physics  data  be  generated 
without  intervention  by  the  user 

The  program  changes  that  accomodate  these  features  would  reduce 
a users  task  to  describing  the  geometry,  source,  and  detector 
of  the  problem. 

The  development  history  of  the  program  file  has  required 
"patching  in"  of  some  physics  models  with  a degradation  in 
efficiency.  The  following  program  modifications  are  recommended 
to  improve  efficiency 

a.  Particle  tracking  be  performed  in  a batch  mode,  i.e., 
perfcm  the  sane  operation  on  a series  of  histories 
before  proceeding  to  the  next  operation 

b.  Precalculate  as  tables  all  primary  to  secondary  transi- 
tion probabilities 

These  modifications  will  also  yield  more  efficient  calculations 
on  vector  arithmetic  computers. 
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APPEI^DIX 


CONTROL  CARDS 


The  users  guide  is  obtained  as  a formatted  listing  using 
the  following  control  cards  on  the  Air  Force  Weapons  Laboratory 
CDC  7600 

TMJAA,T17,STMFZ. 

ACCOUNT  (NA14E,NUlffiER,  GROUP , EXTENSION ) 

ATTACH(  GUIDE,  ID=DYSXTIiJ) 

C0PYSBF( GUIDE, OUTPUT) 


A problem  is  executed  from  the  following  control  cards 


TMJBB.Ti^-OO.STMFZ. 

ACCOUN  T ( NAl'IE,  NUIffiER,  GROUP , EXTENSION ) 
STAGE (MASTER, NT, HD, PRE, VSN=GL97 , ST= ANY) 
COPYP (MASTER, OLDGO ) 

COPYT-  (MASTER,  TAPE9 ) 

RETURN (MASTER) 

REWIND ( TAPE9 ) 

LDSET(MAP=0,FILES=TAPE9) 

OLDGO (PL=1  00000) 

end  of  record 

problem  data 
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